Confined Plasmons in Graphene Microstructures: Experiments and Theory 
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Graphene, a two-dimensional material with a high mobility and a tunable conductivity, is uniquely 
suited for plasmonics. The frequency dispersion of plasmons in bulk graphene has been studied 
both theoretically and experimentally, whereas no theoretical models have been reported and tested 
against experiments for confined plasmon modes in graphene microstructures. In this paper, we 
present measurements as well as analytical and computational models for such confined modes. We 
show that plsmon modes can be described by an eigenvalue equation. We compare the experiments 
with the theory for plasmon modes in arrays of graphene strips and demonstrate a good agreement. 
This comparison reveals the important role played by interaction among the plasmon modes of 
neighboring graphene structures. 



Graphene electronics and optoelectronics have 
emerged as fields of tremendous interest not only 
as improvements to existing technology, but also as 
platforms for completely novel devices. A particularly 
interesting application of graphene is for plasmonic 
devices [TH3], which manipulate charge density waves in 
the two-dimensional atomic sheet. Graphene plasmons 
can have frequencies in the 1-100 terahertz range but 
wavelengths in the micron and sub- micron range [3]-[5] , 
enabling extreme confinement of electromagnetic energy. 
In addition, plasmon frequencies in graphene can be 
tuned through electrostatic [l] or chemical [2] doping, 
making graphene plasmonics a unique platform for 
tunable terahertz sources, detectors, switches, filters, 
interconnects, and sensors. 

The dispersion of plamsons in bulk graphene has been 
obtained analytically [4 , and the experimental results 
have been shown to agree well with the theoretical 
predictions [5]. Plasmons can be confined in patterned 
graphene micro- and nano-structures such that only a 
discrete set of modes can oscillate [1] E]- Such confined 
plasmon modes are of interest for device applications 
since, unlike bulk plasmons, they can couple directly to 
normally incident electromagnetic radiation. No analyti- 
cal techniques have been reported that describe plasmon 
modes in arbitrary graphene structures and model in- 
teractions among plasmon modes of neighboring struc- 
tures. In this paper, we present experimental and theo- 
retical results for confined plasmon modes. We show that 
long wavelength plasmon modes in graphene microstruc- 
tures can be described by an eigenvalue equation. The 
results obtained from the eigenvalue equation match well 
the results obtained using a finite-difference time-domain 
(FDTD) technique. By comparing the measured trans- 
mission spectra of interacting plasmon modes in an array 
of graphene strips with the theoretical results, we show 
that the theoretical models fit the experimental data very 
well. 
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FIG. 1: (a) A cross-section (not to scale) of an array of 
graphene strips with the electric field lines for the lowest plas- 
mon supermode. (b) Optical micrograph of the sample with 
four regions of etched graphene strips (square regions in the 
center) and two reference regions. (c,d) Bright-field optical 
micrographs (lOOX) detailing etched graphene strip arrays af- 
ter the resist was removed. 



A cross-section of the structures considered in this 
work is shown in Figjlja), which shows an array of pat- 
terned graphene strips. The graphene used in our exper- 
iments was grown by chemical vapor deposition (Kevek 
Innovations 1" System) on copper foils and transferred, 
as described by Li et al.[6], to high-resistivity silicon 
wafers (>10 kl^-cm) with ~300 nm of thermally grown 
Si02. Arrays of graphene strips of widths W = 0.75, 
1, 2, and 3 fim were patterned using standard pho- 
tolithography followed by etching in an oxygen plasma 
(see Figjljb)). For all devices, the strip width W and 
the spacing S between the strips were chosen to be equal. 
Graphene strips were doped using HNO3 [7 . The doping 
density was estimated to be in the 4.5-5.0x10^^ cm~^ 
range using the Raman technique described by Das et 
aL[8;. 

Plasmon resonances of arrays of graphene strips were 
measured at room temperature using Fourier-transform 
infrared spectroscopy (FTIR). Figj2] shows the trans- 
mission spectra, of all four strip sizes, for polarizations 
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FIG. 2: (Solid) Measured transmission of radiation polarized 
(a) perpendicular and (b) parallel to graphene strips is plotted 
for four different strip widths W = 0.75, 1, 2, and 3 /xm. For 
all devices, the spacing S between strips is equal to the width. 
A bare Si02/Si substrate is used as reference. (Dots) FDTD 
simulation fits to the measured results. Extracted resonance 
frequencies are 226, 197, 135, and 112 cm~^. 



perpendicular (a) and parallel (b) to the strips. The 
transmission of incident radiation polarized parallel to 
the strips decreases monotonically at long wavelengths, 
showing the expected Drude-like frequency dependence. 
There is no dependence on the strip width. Transmission 
spectra of incident radiation polarized perpendicular to 
the strips show plasmon resonances [TJ [2]. 

The measured transmission spectra can be de- 
scribed qualitatively using a damped harmonic oscillator 
model [H El UHll], 
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Here, r]o is the free-space impedance, a{uj = 0) is the 
DC conductivity of bulk graphene, / is the fill-factor of 
the strips, ngub is the refractive index of the substrate, 
r is the Drude scattering time, Up is the plasmon fre- 
quency, and ^ = (or 1) for incident radiation polar- 
ized parallel (or perpendicular) to the strips. The bulk 
plasmon frequency for small wavevectors is given by the 
expression [4J, 
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where q is the magnitude of the plasmon wavevector and 
is the average dielectric constant surrounding the 
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graphene sheet. Using q = 7r/VF|12|, we find that Eqj2] 
significantly overestimates the plasmon frequencies com- 
pared to the experimental values. For example, using 
Eq{l] to fit the transmission spectra of W = 1 /im arrays 
for parallel polarizations (^ = 0), we find the average 
value of a{uj = 0) and r to be 0.95 mS and 31.5 fs, re- 
spectively. Eqj2]then gives a plasmon frequency of ~245 
cm~^, which is significantly higher than the measured 
plasmon frequency of ~197 cm~^. Such a large error sug- 
gests that better models are needed to understand con- 
fined plasmon modes in patterned graphene structures. 



We first present an analytic technique that captures 
the essential physics of the problem and results in an 
eigenvalue equation for the plasmon modes. Assuming a 
Drude-like conductivity for graphene we start with 
the time-derivative of the equation for the current density 
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Here, Einc is the incident field and is the depolariza- 
tion field that results from the plasmon charge density. 
Only field components in the plane of the graphene sheet 
are included in Eqj3| Using the charge continuity equa- 
tion and the Poisson equation, and ignoring retardation 
effects, the depolarization field can be related to the cur- 
rent density by. 
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The tensor f{f-r') equals [l - 35 (g) s/\s\^] /\s\^, where 
s = f — r' . f{f — f^) is related to the Green's func- 
tion that relates the field to the polarization density 
and can be computed for more complicated geometries 
then considered in this work. We see from Eq|3] that if 
{a{uj = 0)/T)dEd/dt equals —uj^K^ then in the absence of 
£^inc and dissipation the current density will oscillate at 
the frequency ujp. Comparing with Eqj4j it follows that 
the current density associated with the plasmon mode 
satisfies the following eigenvalue equation, 

a(r,a; = 0) j ^2^j^^_ . ^^^^ ^ 2^^^^ 
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The above equation can be solved for the current den- 
sities Km and frequencies Upm associated with the plas- 
mon modes in any graphene structure. The modes satisfy 
the orthogonality condition, J d?fKm{r) • Kp{r) / cr(r, uj = 
0) (X Smp' For the case of bulk plasmons, Eqjslreproduces 
the result in Eqj2] Solving the eigenvalue equation nu- 
merically for the case of a single infinitely long graphene 
strip, we obtain the following result for the frequency of 
the lowest two plasmon modes. 
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The computed current and charge densities for the low- 
est two plasmon modes are shown in Figjsja). Although 
all even plasmon modes (0,2,4...) will couple with nor- 
mally incident radiation, only the lowest plasmon mode 
will couple appreciably. The scaling of the plasmon fre- 
quency with 1 / VW is in perfect agreement with our data. 
For the case W = 1 /im, using Eq|6]and the extracted val- 
ues of (j{uj = 0) and r, ujpo is ~211 cm~^. The eigenvalue 
equation more accurately models the plasma resonance 
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FIG. 3: (a) The computed current densities K(x) (top) and 
charge densities p{x) (bottom) for the lowest two plasmon 
modes of a graphene strip are plotted, (b) The computed 
current densities are plotted for the first 5 supermodes of an 
array consisting of 9 graphene strips {W=S). Locations of 
the strips are indicated by the red horizontal lines. 

than Eqj2] (with q=7r/W)^ but it still overestimates the 
measured resonance frequency by ~6.5%. We next ad- 
dress possible origins of this discrepancy. 

Interactions among plasmon modes in neighboring 
strips can be included by solving Eqjs] using a{f^uj = 0) 
appropriate for an array of graphene strips. Interactions 
lift the degeneracy among strips and result in a band 
of plasmon modes that are the supermodes of the ar- 
ray. The computed current density for the lowest five su- 
permodes of an array containing nine strips is shown in 
Figjsjb) . Only the lowest supermode couples appreciably 
to the normally incident radiation, and it is the frequency 
of this supermode that is measured in our transmission 
experiments. Unfortunately, the matrix eigenvalue equa- 
tion obtained from EqjS] is not sparse, so obtaining so- 
lutions for large arrays is computationally prohibitive. 
Starting from the lowest plasmon mode of a single strip, 
perturbation technique can be used to obtain an expres- 
sion for the frequency ujp{n^N) and the current density 
K{n^N){r) of the nth supermode (n=O...A^ — 1) of an 
A/'-strip array, 

c^2(n,iV)«w2^(^l-2AiCos(^7r^^)^ (7) 
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where Ai is the first nearest neighbor interaction pa- 
rameter. Including second nearest neighbor interactions, 
u;^(0, oo) = cj^Q (1 - 2Ai — 2A2). Ai and A2 are given 
by the expression, 

^ _ /dV/dV^o(r)-/(r-f')-^o(f'-6>(5+iy)f) 
^ ~ / dV/ dV^o(r) • /(r - r') • ^o(P) 

(9) 

Fig|4[a) shows the the calculated values of cj^(0,oo) as 
function of the strip spacing S assuming first and second 
nearest neighbor interactions. The plasmon frequency is 



FIG. 4: (a) The frequencies of the lowest plasmon supermode 
C(;p(0,(X)) of an infinite graphene strip array calculated using 
perturbation theory and FDTD are plotted as a function of 
the ratio S/W . (b) The computed charge densities for the 
lowest plasmon mode in an isolated graphene strip {S/W ^ 
1) and in an array of strips with S/W — 0.25 are plotted. The 
charge densities were computed using FDTD {S/W ^ 1 and 
S/W = 0.25) and the eigenvalue equation (Eq|5]) {S/W > 1). 

reduced as a result of the interactions among neighboring 
strips. For S=W=1 /im, the values of Ai and A2 are .035 
and .009, respectively, resulting in a '^4.5% decrease in 
the value of 6^^(0,00) compared to uj^q. In this case, 
a;p(0, 00) ~ 202 cm~^, which is closer to the measured 
~197 cm~^ than uo^q alone. These results suggest that 
interactions cannot be ignored between nearby graphene 
plasmonic structures. 

The eigenvalue equation does not include retardation 
effects, which could be important in the case of large ar- 
rays, and it also does not account for the discontinuity in 
the field at the oxide/silicon interface (screening by the 
silicon substrate). A technique is needed that incorpo- 
rates these effects, can be used to determine the accuracy 
of Eqjsj and can also compute the measured transmission 
spectra more accurately than Eq{l] For example, Eq{l] 
predicts T^=o(^ = 0) = T^=i(cj = cjp), but the mea- 
sured values in Figj2] differ by ~1.5%. The discrepancy 
arises because the transmission through the gaps in the 
strip array cannot be modeled simply with a fill-factor /, 
especially when the incident radiation is polarized per- 
pendicular to the strips and S < W. 

In the FDTD method. Maxwell's equations are stepped 
in time. In order to model plasmons, we include an aux- 
iliary equation for the graphene current density (Eqjs] 
without the extra time-derivative) and step the equations 
using Yee's leap-frogging algorithm p!4]. This approach 
naturally handles interactions and electromagnetic retar- 
dation. A challenge in the FDTD technique is the range 
of length scales important to the problem. The radia- 
tion frequencies of interest have free-space wavelengths 
extending up to 300 /im, but the corresponding plasmon 
wavelengths are on the order of 1 /im. Furthermore, it 
is important for the modeled graphene thickness to be 
much less than the plasmon wavelength. Therefore, the 
length scales of importance span three orders of magni- 
tude, necessitating a highly non-uniform mesh, with grid 
steps varying from 0.01 — 0.5 /im. The computational do- 
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main spans more than 200x200 /im^, and is surrounded 
by perfectly-matched layer boundaries fT4]. We use the 
values 4eo and 12eo for the THz dielectric constants of 
Si02 and Si, respectively. Plasmonic structures are ex- 
cited at zero angle of incidence with a broadband (0.5-15 
THz) pulse of electromagnetic radiation. The transmis- 
sion spectra for fields polarized parallel and perpendicu- 
lar to the strips are obtained by Fourier-transforming the 
time-domain transmitted pulse. Values of <j(co' = 0) and 
r used in the simulations were iteratively improved until 
the simulated transmission spectra for both polarizations 
optimally fit the measured spectra. 

The FDTD simulation results, shown in Figj2j accu- 
rately fit the measurements. Extracted values of a{uj = 0) 
and r lie in the range 0.91-0.95 mS and 29.5-31.5 fs, re- 
spectively. Using the expression for the graphene con- 
ductivity in Ref.[13 , these values correspond to doping 
densities of 5.0 — 5.2 x 10^^ cm~^, consistent with the den- 
sities determined using the Raman technique. Small vari- 
ations in the parameters across the CVD graphene sam- 
ple are consistent with those measured by terahertz spec- 
troscopy in similar samples[15 . The ability of the FDTD 
technique to quantitatively fit the depth and width of 
the plasmon resonances, while also predicting their cen- 
ter frequencies to an accuracy within one percent, un- 
derscores its usefulness as a tool for modeling graphene 
plasmonic structures. 

The computed x- and ^/-components of the electric field 
for the lowest plasmon supermode are shown in Fig|5] 
for arrays of graphene strips with two different spacings 
(5=2, 0.25 /im and W=2 /im). The locations of the 
graphene strips are indicated by the thin black lines at 
y = 50 /im. The dashed lines indicate the locations of the 
silicon/oxide and oxide/air interfaces. The field is highly 
localized near the graphene sheet, extending a distance 
on the order of the plasmon wavelength. The discon- 
tinuity in the normal {y) component of the field at the 
silicon/oxide interface is also visible. In contrast with the 
5* = 2.0 jam case, when = 0.25 /im, the field in the gaps 
between strips is stronger than the field in the center of 
the strips. This effect helps to reveal the physical ori- 
gins of the interaction between neighboring strips. The 
plasmon charge density that accumulates at the edges of 
an isolated strip generates a depolarization field with 
ujpQ (X \Ed\. In a strip array, this edge charge density is 
partially imaged on the neighboring strips, as depicted in 
Figjlja). This effect increases the depolarization field in 
the gaps between strips but reduces the field within each 
strip. Equivalent ly, the depolarization fields from neigh- 
boring strips are in-phase in the gaps between strips but 
out-of-phase in their centers. Therefore, ^^^(O, oo) < ujpo. 
In contrast, in the highest supermode of the array, the 
current density oscillations in neighboring strips are out 
of phase, so ujp{N, oo) > Upo- 

FDTD simulation can serve to evaluate the perturba- 
tion theory, as shown in Figj4] The two methods are 
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FIG. 5: FDTD simulation results are shown for the x- (a,b) 
and y- (c,d) components of the electric fields for the lowest or- 
der plasmon supermode in two W = 2 /xm arrays of graphene 
strips (5* = 2 /im (a,c), 5'=0.25 /im (b,d)). The dashed lines 
indicate the locations of the silicon/oxide and oxide/air inter- 
faces. 

found to agree when S/W>1^ but not when S/W<^1. 
For example, when S=W in Figjlja), the FDTD calcu- 
lated cjp(0, oo) is lower than ujpo by ~6.5%, in contrast 
with the ~4.5% reduction predicted by the perturbation 
technique. This behavior can be understood by exam- 
ining the plasmon charge density. Fig|4jb) shows the 
FDTD-computed charge density in a strip for the lowest- 
order supermode with S/W^l and S/W=0.25. Also 
shown is the charge density obtained by solving Eq|5]for 
S/W:^l, which is nearly identical to the FDTD result. 
But when 6'/W<Cl, the FDTD calculation reveals that 
the charge density is significantly modified as a result of 
interactions; the charge density increases near the strip 
edges to screen the fields of the neighboring strips. Since 
the perturbation theory assumed that the charge and 
current densities are unmodified from the lowest plamon 
mode of an isolated strip, the results became inaccurate 
when S/W<^1. The good agreement obtained between 
the FDTD method and the analytic model for S/W>1 
suggests that retardation effects do not play a significant 
role in the structures considered in this work. 

To conclude, we presented experimental and theoreti- 
cal results for the confined plasmon modes in graphene 
microstructures. We presented an analytic model which 
captures the essential physics and gives an eigenvalue 
equation for computing plasmon modes. We also pre- 
sented a universally applicable FDTD technique. The 
theoretical models presented show good agreement with 
the measurements, and demonstrate the importance of 
interactions in plasmonic structures. Recently, numer- 
ical and analytical results using other approaches have 
been reported for graphene with periodically modulated 
conductivity p!6HT9] and graphene disk arrays [11] in a 
slightly different context. The present work, to the best 
of our knowledge, is the first time that theoretical and nu- 
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merical models have been presented and tested against 
experiments for confined plasmon modes in graphene mi- 
crostructures. The techniques presented in this paper 
can be used to understand, model, and design complex 
graphene plasmonic structures for applications ranging 
from IR detectors and chemical sensors to plasmonic ra- 
diation sources, oscillators, modulators, and metamate- 
rials. 
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